! This is a 2D advection example using slotted cylinder initial condition and
! double periodic boundary condition for FFSL finite volume scheme.
!
! Li Dong <dongli@lasg.iap.ac.cn>
!
! - 2018-03-24: Initial creation.
! - 2018-03-25: Remove flux splitting (i.e. integer and fractional flux parts),
!               because in 2D this is not very helpful to increase stability,
!               we are still limited by CFL condition.
! - 2019-03-18: Improve readability by add inner and outer operators comments
!               and add flux_x and flux_y arguments to subroutine ffsl.
! - 2024-03-24: Use adv_2d_test_case_mod and put ppm and slope into ffsl_common_mod.
! - 2025-03-23: Add SWIFT splitting to realize fully monotonicity preserving.

program ffsl_adv_2d_case

  use adv_2d_test_case_mod
  use ffsl_common_mod

  implicit none

  real, allocatable, dimension(:,:,:) :: rho   ! Tracer density being advected at cell centers
  real, allocatable, dimension(:,:  ) :: rhox  ! Tracer density due to advective operator along x axis
  real, allocatable, dimension(:,:  ) :: rhoy  ! Tracer density due to advective operator along y axis
  real, allocatable, dimension(:,:  ) :: flxx  ! Flux at cell interfaces along x axis
  real, allocatable, dimension(:,:  ) :: flxy  ! Flux at cell interfaces along y axis
  real, allocatable, dimension(:,:  ) :: flxx0 ! Flux at cell interfaces along x axis
  real, allocatable, dimension(:,:  ) :: flxy0 ! Flux at cell interfaces along y axis
  real, allocatable, dimension(:,:  ) :: divx  ! Divergence component along x axis
  real, allocatable, dimension(:,:  ) :: divy  ! Divergence component along y axis
  integer, parameter :: ns = 2                 ! Stencil width
  integer i, j
  character(30), parameter :: scheme = 'ffsl_2d'

  namelist /params/ nx, ny, nt, dt, flux_type, limiter_type

  call get_command_argument(1, namelist_path)
  inquire(file=namelist_path, exist=is_exist)
  if (is_exist) then
    open(10, file=namelist_path)
    read(10, nml=params)
    close(10)
  end if

  allocate(rho  (1-ns:nx+ns,1-ns:ny+ns,2))
  allocate(rhox (1-ns:nx+ns,1-ns:ny+ns))
  allocate(rhoy (1-ns:nx+ns,1-ns:ny+ns))
  allocate(flxx (1-ns:nx+ns,1-ns:ny+ns))
  allocate(flxy (1-ns:nx+ns,1-ns:ny+ns))
  allocate(flxx0(1-ns:nx+ns,1-ns:ny+ns))
  allocate(flxy0(1-ns:nx+ns,1-ns:ny+ns))
  allocate(divx (1-ns:nx+ns,1-ns:ny+ns))
  allocate(divy (1-ns:nx+ns,1-ns:ny+ns))

  dt = 0.6
  call adv_2d_test_case_init('slotted_cylinder', ns, rho(:,:,old))
  call output(scheme, 0, ns, nx, ny, x, y, rho(:,:,old))

  ! Run integration.
  print *, time_step, sum(rho(1:nx,1:ny,old))
  do while (time_step < nt)
    ! call ffsl_lin96()
    call ffsl_swift()
    call advance_time()
    call output(scheme, time_step, ns, nx, ny, x, y, rho(:,:,old))
    print *, time_step, sum(rho(1:nx,1:ny,old))
  end do

  deallocate(rho, rhox, rhoy, flxx, flxy, flxx0, flxy0, divx, divy)

  call adv_2d_test_case_final()

contains

  subroutine ffsl_lin96()

    ! --------------------------------------------------------------------------
    ! Run inner operators.
    call ffsl(rho(:,:,old), rho(:,:,old), flxx, flxy)
    ! --------------------------------------------------------------------------
    ! Calculate intermediate tracer density due to advective operators.
    call divergence()
    ! Subtract divergence terms from flux to form advective operators.
    do j = 1, ny
      do i = 1, nx + 1
        flxx(i,j) = flxx(i,j) - 0.5 * (divx(i+1,j) * rho(i+1,j,old) + divx(i,j) * rho(i,j,old))
      end do
    end do
    do j = 1, ny + 1
      do i = 1, nx
        flxy(i,j) = flxy(i,j) - 0.5 * (divy(i,j+1) * rho(i,j+1,old) + divy(i,j) * rho(i,j,old))
      end do
    end do
    do j = 1, ny
      do i = 1, nx
        rhox(i,j) = rho(i,j,old) - 0.5 * (flxx(i,j) - flxx(i-1,j))
        rhoy(i,j) = rho(i,j,old) - 0.5 * (flxy(i,j) - flxy(i,j-1))
      end do
    end do
    call apply_bc(ns, nx, ny, rhox)
    call apply_bc(ns, nx, ny, rhoy)
    ! --------------------------------------------------------------------------
    ! Run outer operators.
    call ffsl(rhoy, rhox, flxx, flxy)
    do j = 1, ny
      do i = 1, nx
        rho(i,j,new) = rho(i,j,old) - (flxx(i,j) - flxx(i-1,j)) - (flxy(i,j) - flxy(i,j-1))
      end do
    end do
    call apply_bc(ns, nx, ny, rho(:,:,new))

  end subroutine ffsl_lin96

  subroutine ffsl_swift()

    ! --------------------------------------------------------------------------
    ! Run inner operators.
    call ffsl(rho(:,:,old), rho(:,:,old), flxx0, flxy0)
    ! --------------------------------------------------------------------------
    ! Calculate intermediate tracer density due to advective operators.
    call divergence()
    do j = 1, ny
      do i = 1, nx
        rhox(i,j) = (rho(i,j,old) - (flxx0(i,j) - flxx0(i-1,j))) / (1 - dt * divx(i,j))
        rhoy(i,j) = (rho(i,j,old) - (flxy0(i,j) - flxy0(i,j-1))) / (1 - dt * divy(i,j))
      end do
    end do
    call apply_bc(ns, nx, ny, rhox)
    call apply_bc(ns, nx, ny, rhoy)
    ! --------------------------------------------------------------------------
    ! Run outer operators.
    call ffsl(rhoy, rhox, flxx, flxy)
    ! --------------------------------------------------------------------------
    ! Do SWIFT splitting.
    do j = 1, ny
      do i = 1, nx
        rho(i,j,new) = rho(i,j,old) - 0.5 * (                       &
          (flxx0(i,j) + flxx(i,j)) - (flxx0(i-1,j) + flxx(i-1,j)) + &
          (flxy0(i,j) + flxy(i,j)) - (flxy0(i,j-1) + flxy(i,j-1))   &
        )
      end do
    end do
    call apply_bc(ns, nx, ny, rho(:,:,new))

  end subroutine ffsl_swift

  subroutine ffsl(qx, qy, fx, fy)

    real, intent(in ) :: qx(1-ns:nx+ns,1-ns:ny+ns)
    real, intent(in ) :: qy(1-ns:nx+ns,1-ns:ny+ns)
    real, intent(out) :: fx(1-ns:nx+ns,1-ns:ny+ns)
    real, intent(out) :: fy(1-ns:nx+ns,1-ns:ny+ns)

    real cfl, max_cfl, ql, dq, q6
    real s1, s2, ds1, ds2, ds3
    integer i, j, iu, ju

    max_cfl = -1
    ! Along x axis
    do j = 1, ny
      do i = 1, nx
        cfl = ui(i,j) * dt / dx
        max_cfl = max(cfl, max_cfl)
        iu = merge(i, i+1, cfl > 0)
        select case (flux_type)
        case ('upwind')
          fx(i,j) = cfl * qx(iu,j)
        case ('van_leer')
          dq = slope(qx(iu-1,j), qx(iu,j), qx(iu+1,j))
          fx(i,j) = cfl * (qx(iu,j) + (sign(1.0, cfl) - cfl) * dq * 0.5)
        case ('ppm')
          call ppm(qx(iu-2,j), qx(iu-1,j), qx(iu,j), qx(iu+1,j), qx(iu+2,j), ql, dq, q6)
          if (cfl >= 0) then
            s1 = 1 - abs(cfl)
            s2 = 1
          else if (cfl < 0) then
            s1 = 0
            s2 = abs(cfl)
          end if
          ds1 = s2    - s1
          ds2 = s2**2 - s1**2
          ds3 = s2**3 - s1**3
          fx(i,j) = sign(ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0), cfl)
        end select
      end do
    end do
    call apply_bc(ns, nx, ny, fx)
    ! Along y axis
    do j = 1, ny
      do i = 1, nx
        cfl = vi(i,j) * dt / dy
        max_cfl = max(cfl, max_cfl)
        ju = merge(j, j + 1, cfl > 0)
        select case (flux_type)
        case ('upwind')
          fy(i,j) = cfl * qy(i,ju)
        case ('van_leer')
          dq = slope(qy(i,ju-1), qy(i,ju), qy(i,ju+1))
          fy(i,j) = cfl * (qy(i,ju) + (sign(1.0, cfl) - cfl) * dq * 0.5)
        case ('ppm')
          call ppm(qy(i,ju-2), qy(i,ju-1), qy(i,ju), qy(i,ju+1), qy(i,j+2), ql, dq, q6)
          if (cfl >= 0) then
            s1 = 1 - abs(cfl)
            s2 = 1
          else if (cfl < 0) then
            s1 = 0
            s2 = abs(cfl)
          end if
          ds1 = s2    - s1
          ds2 = s2**2 - s1**2
          ds3 = s2**3 - s1**3
          fy(i,j) = sign(ql * ds1 + 0.5 * dq * ds2 + q6 * (0.5 * ds2 - ds3 / 3.0), cfl)
        end select
      end do
    end do
    call apply_bc(ns, nx, ny, fy)

  end subroutine ffsl

  subroutine divergence()

    integer i, j

    do j = 1, ny
      do i = 1, nx
        divx(i,j) = dt / dx * (ui(i,j) - ui(i-1,j))
        divy(i,j) = dt / dy * (vi(i,j) - vi(i,j-1))
      end do
    end do
    call apply_bc(ns, nx, ny, divx)
    call apply_bc(ns, nx, ny, divy)

  end subroutine divergence

end program ffsl_adv_2d_case
